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ABSTRACT 

The relativistically broadened iron K-a line at 6.4 keV observed in the Seyfert 1 galaxy 
MCG-6-30-15 has provided a probe of the strong-gravity environment near a black 
hole, in particular suggesting that it is rapidly spinning. An important variable in 
such analyses is the geometry of the illuminating source. We present a new technique 
which constrains this geometry based on the spectral line shape, based on a model 
of discrete, point-like flaring regions in the X-ray corona. We apply it to simulated 
reverberation maps and give examples of successful reconstructions of complex coronal 
flaring patterns. For time-averaged spectral lines the problem is highly degenerate, and 
so its inversion more challenging. We quantify this degeneracy and give a measure of 
the spatial accuracy of the method in this case, before checking that it is consistent 
with the existing picture of MCG-6-30-15 by applying it to recent data from XMM- 
Newton. 

Key words: accretion discs - black hole physics - line: profiles - X-rays: general. 



1 INTRODUCTION 

The broad, skewed emission line around 6.4 keV observed in 
the Se yfert 1 galaxy MCG -6-30-15 was first resolved with 
ASCA iTanaka et alllQQsh . A widely accepted model for the 
production of this line is one where hard X-rays are repro- 
cessed by the innermost regions of a cold, thin a ccretion disc 
around a super-m assive black hole llCeorge fc Fabian, 199ll 
iMatt et al.lll99l[ l. X-ray photons emitted from this material 
suffer large gravitational redshifts, special relativistic beam- 
ing and Doppler shifts, combining to give a characteristic 
shape which is in good agreement with the data. In ad- 
dition, a number of alternative models were considered by 
iFabian et al.l (^9^, but a convincing rival candidate was 
not found. 

The extended red wing of this line JWilms et alJl200lt 
iFabian et all 120021: iFabian fc VaughanI l2003h has provided 
the first direct observational probes of the strong gravity en- 
vironment near a black hole and has been used to constrain 
parameters such as the inclination of the accretion disc and 
the b lack hole's spin (jlwasawa et al. 1996; Dabrowski ct al. 
Il997t) . The geometry of the illuminating source of hard X- 
rays, or corona, is a criti cal variable in such calc ulations, as 
illustrated by the work of lDabrowski et all l|l997fl . where the 
authors assumed that the X-ray source was planar, hugging 
a disc which extended down to the innermost stable circu- 
lar orbit in the Kerr metric. Within this framework, they 
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showed that the line profile obtained during a 4 day ASCA 
observation of MCG-6-30-15 in 1994 constrained the spin 
parameter of the black hole to be > 0.94 with about 90% 
confidence. 

The central regions of an accretion disc are, however, 
likely to be subject to a large number of magnetohydrody- 
namical instabilities with the result that a planar geome- 
try for the hard X-ray source is a p oor approximation. An 
altern ative approach was taken by iRevnolds fc BegelmarJ 
1^93). They were able to fit the same data by assuming 
a Schwarzschild black hole where the X-ray source was a 
point-like flare located a few gravitational radii above the 
hole on the symmetry axis of the disc. The bulk of the emit- 
ted radiation would then come from material on plunge tra- 
jectories wi thin the innermost stable orbit radius. Shortly 
afterwards, lYoung et al.l (^9^ showed that a geometry of 
this type predicts a large absorption edge in the spectrum, 
due to the line emitting matter being highly ionised. As no 
such edge has been observed, the data seem to favour the 
Kerr model, although tentatively. 

Of particular interest for the future are time-dependent 
scenarios, where a flash of radiation sweeps over the disc pro- 
ducing line emission in differe nt locations and causes the line 
shape to change with time ("Re ynolds et alJll99gtl . Future 
X-ray satellites such as Constellation- X offer the possibility 
of obs erving these time-depende nt spectra, or reverberation 
maps jYoung fc Reynold jlioool) . However, analyses of such 
data depend on the coronal flaring pattern, and most previ- 
ous work has concentrated on predicting the spectra which 
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result from one or two dominant flaring events in the corona. 
A more likely scenario, however, is one where the line profiles 
in reverberation maps result from an ensemble of overlap- 
ping flares, distributed in time as well as space, and so a 
treatment of more complicated coronal conflgurations such 
as these is vital. Here we aim to develop a technique for 
acquiring knowledge of this pattern from spectral data, con- 
straining the distribution of flares in the corona based on the 
shape of a spectral line. We consider an ensemble of coronal 
flares of arbitrary number, instead of assuming that a single 
flare dominates the emission. 



In general, models are commonly fltted to data via the 
statistic, the best fit parameters corresponding to the 
minimum of a Gaussian distribution. With noisy data this 
can result in over-fitting where many different parameter 
values produce a similarly good fit to the data. In this 
case some form of regularisatio n is often applied, s uch as 
the Maximum Entropy Method (iGuU fc Danielllll97dl . This 
amounts to making some extra basic assumptions about the 
object one is trying to infer, in order to pick out a preferred 
model from a group which are not distinguished by the data 
alone. The exact physics of a magnetically flaring corona is 
complicated, and so rather than trying to treat it directly, 
we can make some progress by making some very general, 
minimal assumptions about its nature. 



Given the idea that individual X-ray flares provide the 
illumination for the production of the iron line, we attempt 
to infer the parameters of a set of discrete, point-like sources, 
which are equally likely to be found at any point over the 
disc. We assume that the number of these sources is Pois- 
son distributed, so that if the data can be satisfied with a 
lower number of flares then it will. Finally, we assume an ex- 
ponentially decaying distribution for the strengths of these 
sources, so that stronger flares are not favoured in the re- 
construction. This set of assumptions corresponds to a tech- 
nique known as 'Massive Inference ' JSkillind f 1998ll and is 
quite general. It can be applied to any situation which can 
be modelled as arising from a number of point source like ob- 
jects, however abstract their physical interpretation might 
be. 



In this paper we begin by introducing the Massive In- 
ference technique and its numerical implementation, in Sec- 
tions Inland 2] In Section 121 we then discuss the physical de- 
tails of our model for the forward problem, that of predicting 
the iron line shape as a function of the coronal geometry and 
other system parameters. After describing the calculation 
with which we compute the line proflle due to a single flare 
we proceed to discuss how these lines can be combined to 
yield a transfer function which maps the parameters of inter- 
est to a composite line profile, in Section [3. II We then show 
some examples of the technique's application to simulated 
data: to reverberation maps in Section |5] and time-averaged 
spectral lines in Section |S| The transfer function for the lat- 
ter case is analysed in Section f6. II before treating data from 
an XMM-Newton observation of MCG-6-30-15 in Section|7| 
Finally, in Section|S|we give our conclusions and discuss the 
application of our technique to possible future data. 



2 MASSIVE INFERENCE PRIOR 

It is standard practice in scientific data analysis to fit the 
parameters of a model to a dataset using a minimisation. 
Often, especially at high levels of noise, a unique minimum 
is not clearly preferred by the data and there is ambiguity 
between different sets of parameters. We can still progress 
in a such a problem though, by making extra assumptions 
about the parameters, which are believed to hold indepen- 
dently of the data. These are encoded in a 'prior' probability 
distribution, which modifies the previous Gaussian distribu- 
tion to one where a clearly defined minimum exists, and so 
a 'best' answer can be found, although it is of course still 
dependent upon the validity of the extra assumptions. 

A well known example of this idea is the 'Maximum En- 
tropy Method' IjGull 19891 . which has been used extensively 
in many areas of physics and astronomy, most notably in 
image processing. Here, the parameters are the fluxes col- 
lected in a set of pixels. By assuming that the fluxes are all 
positive and additive ^, one arrives at the conclusion that 
the image, in the absence of any data, has a prior probability 
distribution given by 



Pr(/i) = cxp(-QrS') 

where 5* is the entropy deflned as 



S{h) 



(hi - nil - hi log (^^) 



(1) 



(2) 



and where TV is the number of pixels in the image, hi (the i"' 
element of the vector h) is the flux in the ith pixel and the 
{mi} represent our 'best guess' at the answer in the absence 
of any data - often a uniform distribution, represents the 
weight we wish to give the Maximum Entropy assumption 
relative to the data. Maximising the entropy essentially cor- 
responds to imposing the least possible prejudice on the im- 
age, and the best image is the one which strikes the optimum 
compromise between satisfying the data, by minimising the 
X^ statistic, and maximising the entropy. 

In this work, instead of using a maximum entropy prior, 
we instead follow the 'Massive Inference' route, which can 
be expected to be particularly appropriate to the case of a 
small number of point-like sources. To provide a rigorous 
framework in which we can combine our prior assumptions 
with line profile data, we use a precise statement of con- 
ditional probability - Bayes Theorem. According to Bayes 
Theorem, a hypothesis h given some data D has a poste- 
rior probability distribution 'Px{h\D) given by the product 
of the likelihood Pr(13|/i) and the prior 'Pi{h) normalised by 
the evidence Pr(Z)) 



Pr(/i|13) 



Pr(/t) Pr(.D|/t) 
Pr(I?) ' 



(3) 



We assume Gaussian random errors, so that the likelihood 
takes the form 



Vr{D\h) oc e 



(4) 



^ Fluxes are necessarily positive, and 'additive' means that the 
flux received across two non— overlapping patches is the sum of 
the individual fluxes. 
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X = 



M 



D-Tl{h)],Hu [D-TZ{h)], 



(5) 



where M is the number of data points in the vector D and 
the covariance matrix H is diagonal because the errors are 
independent. In this work, the hypothesis h is comprised 
of the number, strengths and locations of a group of X-ray 
flares, and the data D are the iron line X-ray fluxes received 
in each energy bin around 6.4 keV. TZ is the transfer function 
which provides a map from the parameter space to the data, 
so that TZ{h) are the (noiseless) 'mock data' corresponding 
to a set of parameters h. The prior probability PT[h{N, x, £)] 
consists of the product of three separate priors for each of 
N (the number of flares), x (an A'^ element vector of flare 
positions) and z (an element vector of flare strengths). 
is assigned either a Poisson distribution 



Pr(7V) 



iV! 



(TV) = a ± 



with mean a, or a wider, geometric distribution 
Pr(iV) = (l-r)r^. 



1 + a 



{N) = a ± ^a{l + a) 



(6) 

(7) 
(8) 



depending on how accurately th e mean number of flares a 
can be estimated (ISkillina^ll999^ . In practice we always use 
the latter, with an additional, exponentially decaying prior 
on a itself. The flare positions are assigned a prior which 
is constant over the allowed locations (in practice covering 
the central regions of the accretion disk) and zero otherwise. 
Their strengths an exponentially decaying distribution 



Pr(z) = [] . 



(9) 



where the decay constants qi are chosen internally based 
on the overall scale of the data, and so do not feature in 
the analysis. This form avoids extremely strong flares being 
favoured in the posterior. 

The full posterior probability distribution can then be 
written 



-a N " -z,79i 



iV! 



Y^[D - 7^(/^)], H^, [D - 7^(/l)],; 



(10) 
(11) 



Finding the value of h which maximises this distribu- 
tion will give the 'best' set of parameters for a given dataset. 
This, however, is only a fraction of the total amount of in- 
formation contained in the answer. In this work we gain a 
knowledge of the full posterior dis tribution by s ampling with 
a simulated annealing algorithm (|Skillinell99fll ) described in 
more detail in Section |1] We can then calculate any moment 
of the distribution, such as the mean and standard devia- 
tion, rather than just finding the best-fitting model. 

As mentioned above, this form of posterior is quite gen- 
eral - in this paper, the 'locations' of the sources are the 
spacetime coordinates of coronal flaring events, but this need 
not be the case. The physical interpretation of the parame- 
ter space is entirely determined by the form of the transfer 



function TZ which provides a map from h to the data (with- 
out noise) 



n-.h^D. 



(12) 



In Section |21 we consider the detailed physics which is in- 
cluded in this object, although the end result of the cal- 
culation is simply a linear mapping from the three objects 
contained in h (ie A'^, x and z) to the data D. 



3 LINE PROFILE CALCULATION 

The 'forward problem' which we attempt to invert in this 
paper is that of predicting the shape of the iron Ka spectral 
line at 6.4 keV. We perform explicit ray-tracing in Kerr grav- 
ity which accounts for all relativistic beaming and Doppler 
shifting effects, using the geodesic equations given in the 
appendix. 

We adopt a 'lamp-post' geometry, where a point source 
orbits above an optically thick, geometrically thin accre- 
tion disc, consisting of particles on circular, time-like pro- 
grade geodesies in the equatorial plane. Its coordinates are 
ts,pa,4's,ha in a cylindrical polar system, where ts and 
(ps are Boyer-Lindquist coordinates and ps and hs are re- 
lated to the Boyer-Lindquist coordinates via = + hi, 
tan^s — hs/ps- Except where explicitly stated otherwise, 
the 'radius' of the source refers to its value of ps, the dis- 
tance from the rotation axis of the black hole. As the source 
is located above the plane, it is not in free-fall, rather, it 
remains above the same point on the disc throughout its 
motion. The precise details of its trajectory are also given 
in the appendix. The disc material orbits with a period given 
by 



r = 4 



(- ^ -3/2N GM 



(13) 



where, in the flrst part of this equation we use conventional 
(S.I.) units (where GM/rc? is dimensionless) , and in the 
second part, f is measured in units of Vg — GM/c? and 
a in units of GM/c. The disc's inner radius is chosen to 
be that of the marginally stable orbit and the disc extends 
outward to 100 rg. We set the black hole's angular momen- 
tum parameter a to 0.998 so that the accretion disc allows 
emission within & r„ llTanaka et al.lll995l:lwilms et al.ll200ll : 
iFabian et al.l2002l:lFabian fc Vaughanl2003ll . In aU the plots 
in this paper, the disc rotates in an anti-clockwise sense. 

We assume the source emits isotropically in its rest 
frame a power-law spectrum with a speciflc luminosity oc 
v~°' erg s~^ Hz"'^. In practice we enforce this by choosing 
the initial directions of the emitted photons such that they 
each pass through the centre of one pixel of a HEALPix'^ 
sphere. We can then compute the incident flux at any point 
on the disc, as measured in the rest fram e of the disc matter 
JLu fc YJ2601I: iMartocchia et al.ll2000t) . to be 



F — 



(l-hZ«d)l+" 



(14) 



^ Hieraxchical Equal Area isoLatitude Pixellisation, see 
|http://www.eso.org/sclence/heaIpbc| 
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where Zsd is the redshift suffered by a photon arriving at 
that point and Ud is the frequency as measured by an ob- 
server co-moving with the disc. To find the flux that leaves 
the disc, we use the foUowing semi-analytical approxima- 
tion for the effect of the photons' incident energy and an- 
gle on the number emitted as iron Ka line radiation, which 
are derived from the Monte Carlo simulations of the de- 
tailed ato mic physics occurring in a cold, non-ionised disc 
jGeoree fc Fabiaft.1991) : 

Nout=g{ei) I N,„{v)f{v)dv (15) 
where 

f{v) = 7.4 X 10"^ +2.5 exp (^-^^^^j (16) 
g{e^) = 6.5 - 5.6 cos + 2.2 cos^ Si (17) 

and where f is measured in keV and ut is the threshold 
for triggering fluorescent emission. Above Vm ~ 30 keV this 
approximation breaks down, but there is a negligible con- 
tribution from this region to the total flux. We obtain the 
following expression for the flux of iron Ka line photons 
leaving the disc 

Fout(^5^^Vd9{e,)j^-^-^-^ /(/.,);.;<i+"'d/., (18) 

We assume that the iron line radiation leaves the disc 
isotropically and is collected in a grid of 2000^ equal area 
pixels (modelling those of a CCD detector) in an asymp- 
totically flat region of spacetime (in practice 1000 rg), at a 
location inclined at 30° to the disc, in keeping with appropri- 
ate parameters for MCG-6-30-15. To calculate the change 
in intensity of radiation on its journey from the disc to 
the observer we use the Liouville invariant (e.g. Misner et 
al. 1973) 

— constant, (19) 

where the change in frequency is computed for each pixel in 
the detector by integrating photon trajectories backwards in 
time until they hit the disc or are lost into the black hole 
jLaoilll99ll: iDabrowski et al.lll997^ . The power collected in 
one pixel in the flnite frequency range is then 

Wa. cc J Fout{rA) (20) 

where we have made the position dependence of the emis- 
sivity over the disc explicit and assumed isotropic emission 
from the disc. We integrate over the entire solid angle sub- 
tended by the disc at the observer by summing the contri- 
bution from each pixel. The contribution to the total flux of 
each photon is then binned in both energy and arrival time, 
to build up a time-dependent spectrum. 

3.1 Transfer functions 

By following the procedure outlined in the previous section, 
we can generate the time-dependent spectrum produced by 
a single flaring event at a given location in spacetime. We 
can now obtain the transfer function TZ which maps the co- 
ordinates of an arbitrary number of flaring events to the 
data D. If the corona itself is optically thin (so that it has 



a negligible effect on the iron line emission), then we can 
simply superpose the spectrum coming from each individ- 
ual flare, accounting for the time delays due to both the 
different source positions and the different time at which 
each flare becomes active. In this way, the transfer function 
becomes linear, although this is not a requirement for the 
inversion algorithm to work. 

In this work we concentrate on the radial and azimuthal 
locations of the sources rather than their heights. Rather, we 
fix the height at a value which corresponds to a characteristic 
scale height of the corona which in practice is set to 2 or 
5 r-g , although it would be simple to extend these models to 
include the heights of the fiares as well. The time lag arising 
from a difference in the ts coordinate of the source is handled 
trivially by simply shifting the total spectrum obtained from 
a single fiare along the time axis by an appropriate amount. 
The delays coming from different values of ps and (j)s are 
handled by finding the earliest time at which photons from 
a particular source arrive at the observer for each source 
position. The entire spectrum is then shifted by this amount 
in the full transfer function. Figure [l(a)| shows the resulting 
spectrum when the transfer function acts on a set of four 
fiares, the coordinates of which are given in the caption. 
These spectra s how the same qualitative features as those 
in the work of iRevnolds et a l. (1999), who gave the first 
examples of such time-dependent spectra. 

While the aim of this paper is to present an inver- 
sion method for reverberation maps, data of this nature 
is not yet available. Rather, the best example of a broad 
iron line to date, seen in MCG-6-30-15, requires integration 
times so long that only basic variability analyses are viable 
jlwa^awa et a.] J 1199^: IVaughan EdelsonlEoOll- IShih et a.lJ 
l2002h . In this case to get a mean error of ~ 15%, ob- 
servations of hundreds of kiloseconds have been necessary 
(J'abian ct al. 2002). A recent estimate of the mass of MCG- 
6-30- 15 gives an upper limit of ~ IO'^Mq llMorales fc Fabiaiil 
l2002h . corresponding to a characteristic time-scale {GM/c?) 
of 49 seconds. For radii in the range ~ 2 — 15, Equa- 
tion ll3l implies orbit time-scales of less than one hour, which 
means that observations made with current telescopes such 
as XMM-Newton and Chandra span many orbital periods. 
This clearly calls for a modification of the technique as 
stated above, which we can accomplish by integrating over 
both ts and (j)s to yield a time-averaged spectrum for each 
radius of a fiare, modelled as the ring it sweeps out over 
time, assuming it persists for an entire orbit. In practice, we 
can average the emissivity Foutir, 0) over ij> to obtain an ax- 
isymmetric profile and use a single time bin extending over 
the full duration of fiux detection. In this case we can plot 
the entire transfer function, as shown in Figure [l(b)| 

In reality, of course, we do not expect a perfectly ax- 
isymmetric system, due to the turbulent and non-linear na- 
ture of the disc and corona. To begin to probe these effects 
one can imagine an intermediate scenario, where a time- 
averaged spectrum is obtained over an interval much shorter 
than the orbit time for a fiare, or one where fiaring events 
occur in 'fiashes' which persist for times much less than the 
orbital period. In this case, we lose all the temporal informa- 
tion in the system but we can retain the azimuthal positions 
of the sources. The transfer function then consists of a time- 
averaged spectrum for each value of the pair Ps,4's taken 
by a fiare. In Figure [l(c)| we show spectra for a source at a 
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Figure 1. (a) Time-dependent spectrum produced by four flares: three simultaneous events at {ps,(j>B) = 
(19rg,0°), (lOrg, 100°), (19r9,240°) followed by a repeat of the latter flare 4tg(= GM/c^) afterwards, (b) The fuD '1-d' trans- 
fer function for time-averaged line profiles, generated from azimuthally averaged emissivity distributions corresponding to the 
ring-shaped flares swept out by orbiting flares, (c) Time-averaged spectra for a flare at ps = 1.5 rg for each of the 36 azimuthal bins in 
the '2-d' transfer function, (d) Same as (c) but for ps = lOrg. 



small radius ps = 2, for all values of Figure [T(d) [ shows the 
equivalent plot for a larger value of ps- Averaging the spec- 
tra from these graphs over (j)s would produce cross-sections 
of Figure [l(b)| at the corresponding radii. 



4 COMPUTATIONAL METHOD 

The eventual goal of this work is to construct a framework 
for determining patterns of coronal flaring activity from iron 
Kq line spectra. We can assess the validity of conclusions 
drawn within this framework by testing the technique with 
data generated from a known distribution of flares. We op- 
erate on this distribution (known as the 'truth') with the 
transfer function TZ to yield the 'mock data' TZ (jij . To com- 
plete the process we add a known amount of Gaussian noise 
0", to obtain the simulated data D = TZ (ji^ + a. 

Given these simulated data we can construct an infer- 
ence h via the Massive Inference method, which can then 



1. Choose h ^ 2. Act with Tl 

compare -|- noise 

Mass. Inf. 

4. Obtain inference {h) -< 3. Add noise 

Figure 2. A schematic representation of the process by which 
we tost the technique 



be compared with the truth h. A graphical summary of this 
process is shown in Figure |5| 

Sampling of the posterior is achieved with a Markov 
Chain Monte Carlo algorithm, using a simulated annealing 
cooling schedule. Each consecutive sample is chosen accord- 
ing to a 'temperature' parameter, which is increasingly re- 
strictive in terms of the changes it allows as the system cools. 
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the transfer function on the mean of the sampled posterior. As the algorithm progresses, it can be seen to converge toward the true 
posterior, arriving at a of 48.8/50. (b) A snapshot of some of the final samples. The true posterior is shown as the contour plot and 
the small circles represent the samples, showing that they are concentrated in regions of maximum posterior probability. The diagonal 
reflection symmetry is due to the flares being indistinguishable. 



The samples are initially spread over a large region of pa- 
rameter space, but as the temperature drops, they become 
more concentrated around regions of high posterior proba- 
bility. The algorithm can store an ensemble of several flare 
configurations, each evolving independently, and the system 
is cooled as slowly as computation time allows. In this way, 
we maximise the number of samples taken and minimise the 
chance that the samples are trapped within local minima. 

There is no formal convergence test for the Massive In- 
ference method. Rather, one can vary both the rate at which 
the algorithm anneals and the size of the ensemble of simul- 
taneously evolving parameter sets, and check that they all 
converge to a single distribution. A representation of con- 
vergence for a simple example involving only two flares is 
shown in Figure |3] The flares are located at radii of 5 and 
13 rg. Figure |3(a)| shows the misfit between the mock and 
simulated data as the algorithm progresses from its initial, 
hot state to the final cool state. The decrease in tempera- 
ture is visible as the tapering of the curve's overall envelope 
of variation, while the convergence is reflected in the ap- 
proach of (per degree of freedom) toward the final value 
of 48.8/50. In Figure |3(b)| we plot a projection of the sam- 
pled posterior into a plane showing the radial coordinate of 
each flare, together with a contour plot of the true posterior, 
showing that the samples are indeed concentrated around its 
peak. 



5 TIME-DEPENDENT INVERSIONS 

We now present some examples of the application of our 
technique to simulated reverberation maps, obtained by 
adding noise to spectra such as that shown in Figure [l(a)| 
We quote noise levels as mean signal-to-noise ratios (SNR), 
the ratio of the noise to the average over all pixels in which 
the flux is non-zero. For a particular example of simulated 



data with a mean SNR of 1, which we take as our starting 
point for illustrating the technique, the inferred number of 
flares is given in Figure |4(d)| where a strong peak is visible 
at the true number of 7. The locations of the flare ensem- 
ble above the accretion disc (their Ps,0s coordinates) are 
shown in Figure |4( a) | where each sample of the posterior is 
plotted as a set of grey circles, where A'^ is the number 
of flares corresponding to that particular sample. In reality, 
the inner 20 Vg of the disc are covered by a 20 radii x 18 az- 
imuths polar grid of flare positions, but to help visualise the 
results, we have displaced each grey circle in a random direc- 
tion from the centre of its grid cell by a radius obtained by 
sampling from a Gaussian distribution with a standard devi- 
ation equal to the radius of the thicker, black circles, whose 
centres mark the location of each flare in the true conflgu- 
ration. The radius of each grey circle is proportional to the 
strength of the flare for that sample and is normalised such 
that the strongest flare's radius matches that of the thicker 
circles marking the truth. 

In addition, we plot two circles centred on the origin, 
which mark the inner radius of the accretion disc and the 
radius ps = 20 Vg which corresponds to the outer limit of 
the polar grid, but not the accretion disc itself {vmax = 
100 Tg). The samples are clustered around the true positions 
in all cases but one, showing that the system copes well with 
a problem of this complexity, at this level of noise in the 
data. Figure |4 ( b ) | shows the inferred distribution of times at 
which each flaring event occurred, as a projection onto the 
t-x plane (the observer is located above the negative y axis, 
nearest the bottom of the plot). For the flare closest to the 
black hole whose x-y position was least accurately inferred, 
the samples' time coordinates are centred around the true 
value of 6 ig. It is very difficult to pick out 7 distinct spectra 
by eye in the simulated data shown in Figure [4(c)| and yet 
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-20 -15 -10 -5 5 10 15 20 
Distance (GM/c^) 

(a) 




"0 5 10 15 20 25 30 



Time (GM/c ) Number of Flares 

(c) (d) 

Figure 4. (a) Samples from the posterior probability distribution for a configuration of 7 flares, projected into the equatorial plane. 
Only 6 flares are visible because two are located above the same point in the disc, but are separated in time. The largest circle marks 
the edge of the polar grid of allowed positions, not the edge of the accretion disc, (b) The same samples, but now projected into the 
t-x plane. The observer is located above the negative y axis, nearest the bottom of (a), (c) The simulated noisy data from which the 
inference was made, with a mean SNR = 1. (d) A histogram showing the inferred number of flares, A'^. 



the algorithm is able to do so, and give the locations of the 
flares with a good degree of accuracy. 

We show two more examples in Figures |3 and |S| one for 
a configuration of 11 flares and a mean SNR = 2, and one 
for poorer quality data (SNR = 0.5) with fewer flares in the 
truth. Figure |5(a)| shows a 3-d representation of the locations 
in spacetime of the true and inferred distributions of flaring 
events. The vertical axis measures time, not the height of 
the flares, which is fixed at 5rg. The light grey samples on 
the disc and at the side of the plot correspond to the x-y 
and t-x projections like those shown in Figures [4 (a)| and |4(b)| 
for the case of seven flares. The truth is plotted here as a 



set of wire-frame cubes with a side length equal to twice the 
standard deviation of the Gaussian distribution from which 
the amount of radial displacement of each dark grey sample 
was taken. 

The noisy data and a histogram of the inferred number 
of flares are shown in Figures |5(b)| and |5(c)| It is clear that 
at least 11 flares were chosen each time in order to flt the 
data and higher numbers become increasingly unlikely, re- 
flecting the prior probability distribution on TV. Information 
from Figures [5(a) | and |5(c)| can be combined by plotting the 
latter in colour. Instead of grey circles, we can use a different 
colour for each value of A*' in the histogram. We avoid cer- 
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Figure 5. (a) Samples from a posterior distribution corresponding to a truth consisting of 11 flares, shown as wire-frame cubes. The 
time at which each flare became active is plotted along the vertical axis. The height of each flare is fixed at hs = 5rg. The lighter circles 



plotted on the equatorial plane and the vertical plane at the back of the plot correspond to the two projections of Figures |4(a) | and |4 (b) | 
for the case of 7 flares, (b) Simulated reverberation map from which this inference was made, at a mean SNR = 2. (c) A histogram 
showing the distribution of values of N, the number of flares, chosen by the algorithm during sampling. 



tain colours being obscured by others plotted over the top by 
randomising the order in which they are drawn, thus giving 
a visual impression of the histogram, which can be useful in 
revealing spatial bias of different values of A^. Colour ver- 
sions of all relevant plots in this paper are available online'^ . 
Figure |S| shows an inversion performed for data at a higher 
noise level (mean SNR — 0.5). In both 3-d plots, the thick 



http: / /www. mrao.cam.ac.uk/~rg200 



arrow drawn in the plane of the disc is directed towards the 
observer. 

As the level of noise is increased, the quality of the in- 
ferences decreases - flares begin to be missed or show up in 
incorrect places. Rather than perform a detailed analysis of 
the limits of performance of the above time-dependent setup, 
we now consider those of time-averaged spectra, where we 
focus on just the pa and 4>s coordinates of the flares. The 
reason for this two-fold. Firstly, the behaviour of the system 
for demanding noise levels in time-averaged data is a see- 
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Figure 6. (a) Samples from the posterior distribution for a truth with 3 flares, made at a higher level of noise (mean SNR = 0.5) in the 
simulated data, shown in (b). For this inference, the histogram for A'^ (not shown) shows a single dominant peak at the correct number 
of 3 flares. 



nario more likely to be of interest in the short to medium 
term, given the quality of current data. Secondly, obtaining 
the locations of a set of flares in time is a relatively easy 
problem compared to that of revealing their spatial loca- 
tions, simply because the structure induced in the data for 
a single flare at different times is far more distinctive than 
at different positions. It was described in Section l3.ll that 
the former is achieved by simply translating the spectrum 
from one flare along the time axis in the data-space by an 
amount corresponding to its time-delay, whereas the spec- 
tral differences that come from changes in spatial location 
are far more subtle. It is therefore clearer to focus on the 
time-averaged scenarios described in the next section. 



6 TIME-AVERAGED INVERSIONS 

We can recast the above system in a form which does not 
depend on time by summing the spectra over their entire du- 
ration, and removing all reference to time from the transfer 
function. In this way, we obtain a 2-dimensional problem, 
where we aim to recover the ps and (f>a coordinates of each 
source, which take values in the range 1.5 - 15 rg and - 
360° over a grid of 40 radii x 36 azimuths. In doing this, we 
have reduced the number of degrees of freedom (d.o.f.) in 
the data by almost two orders of magnitude, and as a result 
we readily find noise levels at which the system struggles to 
produce correct inferences. 

Figure |7(a)| shows the set of samples of the posterior 
probability distribution for an inference based on the data 
of Figure |7(b)| The true distribution contains three flares at 
the locations shown by the thick black circles in Figure |7(a)| 
In addition to being spread out over many cells of the polar 
grid of different pa, (j)s values, the inference is also clustered 
around a region above the rear, receding part of the disc, in 
which there is no flare in the truth. The mock data do pro- 
vide a satisfactory fit to the simulated data, however, which 



implies that, at this level of noise, we are probing degener- 
ate regions of the transfer function. We note further that the 
samples of N shown in Figure |7(c)| peak very strongly at the 
true number of three, so that the extra cluster of samples 
in Figure |7(a)| does not imply the presence of a fourth flare, 
rather an additional possibility for the location of one of the 
three. 

6.1 Linear-dependence in the transfer function 

We can quantify the linear-dependence present in the trans- 
fer function by performing a singular value decomposition, 
yielding the singular values marked by the squares in Fig- 
ure Is^aj] The gradient of this line determines the decrease 
in noise that is necessary, on average, in order to detect 
the presence of extra flares as they are added into the true 
distribution. The 'flare-space' singular vectors (in the space 
of all possible flare configurations) always form a linearly- 
independent basis, with which we can represent any flare 
pattern. However, their relative magnitudes are given by 
the singular values, which means that for a given finite noise 
level, a certain (large) fraction are effectively zero, remov- 
ing this ability to span flare-space and represent an arbi- 
trary pattern. Degenerate regions arise when the projection 
of some of a particular flare pattern lies in these zero direc- 
tions and fails to result in detectable features in the data. 

A useful measure of the overall effect of this is shown 
in Figure |8(b)| where we have summed each singular vector, 
weighted by the singular values, and projected the result 
onto the polar grid covering the inner regions of the accre- 
tions disc. A grey circle is drawn at the centre of each grid 
cell, whose shading gives the overlap of the average singular 
vector with the cell. The radii of the circles increase linearly 
with radius to provide a suitable covering of the disc, for 
ease of visualisation only. The plot reveals those cells which, 
on average, produce the same type of structure in the data 
and which require higher levels of noise to distinguish be- 
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Figure 7. (a) Posterior samples for a distribution of 3 flares, inferred from time-averaged data, hence only an equatorial plane projection 
is available. The true flare positions are marked with thick circles and the large circle at 15 Vg marks the outer limit of allowed flare 
positions, (b) The time-averaged spectral line used in this reconstruction, shown as the points with error-bars. The mean signal-to-noise 
ratio was 7, corresponding to an average noise level of 15%. The solid line shows the mock data from the inference, (c) A histogram 
for the value of TV corresponding to each set of samples shown in (a). Although the samples are clustered around four distinct regions, 
the algorithm rarely picked sets of four flares, instead preferring the correct number of three. 



tween them. This explains the patterns seen in the inferred 
flare locations shown in Fi gure |7(a)| as the algorithm ex- 
plores the contours of Figure |8(b)| with approximately equal 
frequency, in its search to fit the data. This plot also re- 
veals the effective spatial resolution of the technique, as a 
function of position on the disc. Two adjacent flares placed 
in cells corresponding to a high gradient in Figure |8(b)| will 
be easier to resolve than two placed in a region where the 
gradient is small. 

As described in Section [3. II we can reduce the level of 
complexity in the transfer function still further by calculat- 



ing the time-averaged spectra which arise from axisymmetric 
distributions of emissivity due to 'ring-shaped' flares, mod- 
elling time-averaged point-like flares which persist over an 
entire orbit or more. We generate these axisymmetric emis- 
sivity profiles by averaging the illumination pattern due to 
a single flare over all values of (ps , resulting in a 1-d problem 
where we attempt to recover just the radius ps of the ring. 

Performing this procedure results in a transfer function 
whose singular values are given by the circles in Figure |8(a)| 
They drop much more quickly than for the 2-d version, giv- 
ing a difference of one order of magnitude for just 5 fiares. 
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Figure 8. (a) Singular values of the time-averaged transfer functions, for the 2-d system (inferring ps and 0a) shown as squares and the 
1-d system (inferring just ps) shown as circles. Losing azimuthal information results in a dramatic reduction in the transfer function's 
linear independence, (b) The average of the flare-space singular vectors for the 2-d system, showing degenerate regions above the accretion 
disc, within which a flare produces very similar features in the data. On the approach (left-hand) side of the disc, beaming produces 
sharper structure in the data which leads to smaller degenerate regions and a greater spatial accuracy, (c) Inferred (ring-shaped) flare 
radii and strengths for 3 flares at different levels of noise, plotted as the mean flare strength in each of the 40 radial bins. At SNR = 10, 
two flaring regions are revealed, but a further factor 10 reduction in the noise is necessary in order to successfully resolve the 3 flares. 



This trend is reflected in the picture of the 1-d transfer func- 
tion shown in Figure [l(b)| where the line profile shows only 
modest variation as a function of radius. To illustrate the 
difflcuhy in performing inference calculations with this sys- 
tem, we show the results of many such calculations for a 
wide range of noise, in Figure |8(c)| where we plot the mean 
strength of the samples in each radial bin. At high noise lev- 
els (SNR ~ 1) none of the three flares present in the truth 
was located in the inference. For noise levels comparable to 
the 2-d calculation of Figure |7| the innermost flare was de- 
tected, but it required a further drop of about a factor of 10 
in noise before all three flares could be detected. 



7 APPLICATION TO MCG~6-30-15 

We now apply the technique to the results of a 325 ks obser- 
vation o f MCG-6-30-15 mad e with XMM-Newton, described 
fuUy in iFabian et al.l ijioO^) ■ We take the spectral line to 
be the difference between their Model 4 for the underlying 
continuum, and the data. The average noise level is approx- 
imately 15%, which corresponds to a mean SNR ~ 7. Given 
the results of Figure |8(c)| it seems highly doubtful that any 
meaningful inferences can be made within the 1-d axisym- 
metric, time-averaged framework at the level of noise. We 
choose instead to apply the 2-d system, but with the fol- 
lowing large caveat. Recent estimates of the mass of MCG- 
6-30-15 give values close to 10^ Mq giving a characteristic 
time-scale of 49 seconds, which, as discussed in Section f3. II 
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Figure 9. (a) Inferred distribution of samples for the XMM-Newton data shown as the error-bars in (b). (b) The solid line shows the 
mock data obtained by operating on the inferred flare distribution with the transfer function, together with the data themselves. 



corresponds to an orbital period of less than one hour even 
for a relatively slow flare, orbiting at — 15 r^. 

It is clear from these time-scales that we cannot hope 
to attribute structure in the line profile to particular non- 
axisymmetric features in the results of any inference cal- 
culation, although such results may be useful in giving an 
indication of the general radii at which flares could be ac- 
tive. In addition, we make no attempt in this work to model 
the form of the underlying continuum spectrum, and so are 
completely reliant on such a model to define a given spec- 
tral line shape. Models for the continuum, while providing 
a good overall fit, still allow considerable variation of the 
slope in the vicinity of the iron line, thus altering the shape. 
Detailed predictions of flaring activity based on the precise 
details of this line profile can therefore be easily changed 
with a different continuum model, and should not be taken 
as clear evidence of a particular flare distribution. This pa- 
per is intended to demonstrate the type of calculation which 
will become possible with the advent of high quality, time- 
dependent spectral data. The application of our techniques 
to this XMM-Newton observation is a check, to make sure 
we obtain a plausible picture from current data which is 
consistent with previous studies of this object. 

Figure |9(a)| shows the results of applying the 2-d sys- 
tem to these data. The posterior samples are concentrated 
in two regions, and the histogram for the number of flares 
(not shown) has a very strong peak at = 2. The solid 
line in Figure |9(b)| shows the mock data from this inferred 
flare pattern, which provides a good fit with the data below 
~ 7 keV. The structure abo ve ~ 7 keV is due to iron K/3 
emission llFabian et alJ2002!l , which we do not include in our 
transfer function (although it would be simple to do so). For 
this reason, we ignore the points above 7 keV when calcu- 
lating the value of for this fit, which is 0.98 per degree 
of freedom. That the flaring regions are contained within a 
few gravitational radii of the central black hole, which leads 
to a very steep emissivity profile in this region, is consistent 



with the existing picture of this object based on previous 
studies. 

Whether the precise azimuthal distribution of flares has 
a direct physical signiflcance, given the time-scales involved, 
is doubtful. However, it is important to assess the unique- 
ness of this prediction - to what extent would the result 
change with a different realisation of the noise in the data? 
To address this issue, we can perform a further simulation 
based on the results of Figure |9(a)| as follows. By taking 
a truth consisting of the mean flare-strength in each grid 
cell, we can generate mock data to which we add an amount 
of noise corresponding to that in the original data in Fig- 
ure |9(b)| We can then perform repeated applications of the 
inversion technique for different random seeds in our noise 
generator. Figure IIUI shows the results of two such calcu- 
lations. The thick circles show the 'truth' consisting of the 
results of Figure |9(a)| and the grey circles show the samples 
of the posterior probability. In Figure [lO(a)| we obtain very 
similar results to the original inference, but in Figure [lO(b)| 
we see some differences. Again, we can interpret this result 
within the singular vectors of Figure |8(b)| - the 'erroneous' 
samples are continuously connected along regions of high 
degeneracy. For this reason, we should interpret the results 
of Figure |9(a)| as indicating the presence of flares anywhere 
within the two regions in Figure |10(b)| Figure [To(cyi shows 
the noisy data for the two different realisations of noise and 
the solid and dashed lines are the mock data correspond- 
ing to Figures |10(a)| and |10(b)| respectively. Another way 
of seeing how similar the two different flaring patterns are, 
when considered in the context of the singular vectors, is 
to evaluate the components of the two patterns on the sin- 
gular vector basis. The solid line in Fig ure |10(d)| gives the 
projection of the 'real' results of Figure |9(a)| onto the sin- 
gular vectors, and the circles and dotted line correspond to 
those for Figures [lO(a)| and |10(b)| respectively. It is through 
projections like these that we can assess the true similar- 
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Figure 10. (a) Samples from the results of applying the 2-d system to noisy data generated from a truth consisting of the results of 
Figurc |9(a)| using a noise level equivalent to that in the XMM-Newton data, (b) Same as (a) except for a different random seed used to 
generate the data, (c) The data themselves. The squares show the data used to obtain (a) and the circles (b), with the error-bars omitted 
for clarity. The solid line shows the mock data from the results in (a) and the dashed line for (b). (c) The projection of the truth (ie 
the results of Figure |9(a)^ onto the first twenty singular vectors, shown as the solid line. The remaining coefficients are negligible. The 
circles show a similar projection for the inference in (a) and the dotted line shows the projection of the inference in (b). 



ity of structure in the data that different patterns of flaring 
activity can produce. 



8 CONCLUSION 

We have presented a new way of analysing line profile 
data which constrains the number, positions and strengths 
of hard X-ray sources in an accretion disc's corona. The 
method makes use of the 'Massive Inference' prior probabil- 
ity distribution which is based on an assumption of discrete, 
point-like sources. 

By applying it to simulated noisy reverberation maps, 
we found that the technique performs well, making the most 
of the rich information content of such data. The algorithm 



was able to disentangle the complicated signatures of combi- 
nations of flare events, rather than concentrating on scenar- 
ios where it can be assumed that one or two dominant flares 
are responsible for the ob served data iRevnolds et alJll999l: 
lYoung fc Revnoldjl200(]ll . Indeed, if there is only a single 
large flare present, such a conclusion will be reached by our 
method. 

Reverberation mapping of relativistic spectral lines is 
still some way off, and so we modified our system to treat 
time-averaged spectra at similar levels of noise to those ob- 
tained by current instruments such as XMM-Newton and 
Chandra. In doing so we lost all the temporal information 
in the transfer function and reduced the problem to that of 
inferring the set of radial and azimuthal coordinates of the 
flares. The relative information loss in the data, however. 
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was greater than in the parameter space, and the resulting 
inference problem was more challenging. At realistic levels 
of noise, we found some ambiguity in the inferred spatial lo- 
cations of the flare, but less in the inferred number present 
in the ensemble. This ambiguity was quantified using a sin- 
gular value decomposition, which yielded a 'map' of the disc 
showing the relative ease with which different regions of the 
disc could be differentiated, and gave a limit on the spatial 
resolution of the technique. 

We reduced the problem still further, considering data 
collected on time-scales greater than the orbital period of 
each flare. In this case we averaged the illumination pattern 
of a flare over the azimuthal coordinate (j)s to obtain a 1-d 
transfer function based on 'ring-like' sources, and attempted 
to constrain the radii of these rings. This problem proved to 
be very difficult, requiring very low levels of noise to ob- 
tain useful results. This very high level of degeneracy in the 
transfer function is due to the line profile being affected only 
weakly by the position of the source, the dominant contri- 
bution coming from the relativistic effects occurring during 
the photons' subsequent journey from disc to observer. 

As a consistency check, we applied the technique to 
a recent observation of the iron Ka line in MCG-6-30-15 
made by XMM-Newton over 325 ks. This long integration 
time means that, given the mass of the central black hole 
of ~ lO^M©, it is not possible to resolve any structure in 
the data due to non-axisymmetric illumination of the disc, 
and that the f-d ring-based transfer function is most appro- 
priate. However, given the difficulties present in the highly 
degenerate ring-based system described above, we applied 
the 2-d time-averaged algorithm to this data, allowing some 
azimuthal structure in the corona. We found two regions of 
flaring activity, one at a radius of ~ 2 on the near side of 
the disc which is responsible for the bulk of the highly red- 
shifted emission below 5 keV, and one on the approaching 
side at ~ 6rg which is more beamed and produces the tall 
peak close to 6.4 keV. 

Further simulations for different realisations of random 
noise illustrated explicitly that this result must be inter- 
preted with the singular vector map of Figure p(b)| in mind, 
spreading out the regions occupied by the samples to cover 
those in Figure [To(b)| This approach showed that, within the 
assumptions underlying the method (mainly that different 
locations of X-ray flares in the strong-gravity region of the 
central black hole account for structure in the line profile), 
the results of Figure |9(a)| are quite robust. Physically, how- 
ever, we should not expect to be able to pick out any partic- 
ular azimuthal region based on a 325 ks observation. In fact, 
by including a relatively weak, narrow component in the line 
in order to account for the tall peak near 6.4 keV, it is pos- 
sible to achieve a very good fit with the da ta based on an 
azimuthally averaged illumination pattern jMiniutti et alJ 
[2003). In addition, the uncertainty that comes from defining 
a spectral line shape relative to a given model for the un- 
derlying continuum radiation, means that the above results 
(other than the loose condition that the fiares are contained 
within ~ 15 Vg) can be changed very easily by making differ- 
ent assumptions about the continuum. For these reasons, we 
do not claim that the results of Figure |9(a)| are a meaningful 
prediction for MCG-6-30-15, rather, it is reassuring to ob- 
tain a pattern contained within a few gravitational radii of 
the central black hole, and gives an indication of the sort of 



results which may become viable with higher quality data. It 
is with the advent of high-quality, time-dependent spectral 
data, like that promised by missions such as Constellation- X, 
that inversion techniques such as those described here offer 
the possibility of mapping the hard X-ray emission within a 
few gravitational radii of a black hole. 
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APPENDIX A: GEODESIC MOTION 



The gravitational calculation s in this paper were performed using geometric algebra, within a new formulation of gravity 
theory recently developed bv iLasenbv et all (ll99Sf l. In this appendix we translate some of the results of these calculations 
into conventional notation, but mention only in passing the motivation for casting the problem in this form. For a detailed 
accou nt of the tools and ideas which lead to these results, see Lascnbv ct al. (1998), Hcstoncs (1966) and Hcstcncs fc Sobczv^ 



We begin with an orthonormal tetrad with components on the Boyer-Lindquist coordinate frame given by 



2 I 2 
t = ^ Ot + 



dr 



■de 



a sin 6* 



(Al) 



p%/A pVA ' p " p"" P ' ps\n0 

where 

= + cos^ e, A = + - 2r. (A2) 
We can then parametrise a photon's momentum with the direction-cosines a and /3 and the energy parameter $ as follows 
p = $ (t + cos a cos 13 f — sin a cos [5 9 ~\- sin /3 0) . (A3) 

This form is based on the (principal) null vector t + f , on which we perform two orthogonal rotations by the angles a and /3. In 
this way, the vanishing norm of the vector is built into its parametrisation, removing a degree of freedom from the calculation 
from t he outset. In addition, we obtain a simple picture of the momentum vector in terms of rotors (see e.g. iLasenbv et all 
il998h \ This form of p yields the following first derivatives of the coordinates x^: 



t = 



(asin6lsin/3VA + r + a 



<i>Va 



cos a cos /3 6 



$ . $ 

— sin a cos /3 </> = — 
P P 



pVA ' ' p 

Substituting p into the geodesic equation results in the following expressions for d and (3 



a sin B 

-F= + 

/A sm9 



$ = 



cos/3 



2p3 
p^ cos /? 



a sin 26 sin a + ( — 4ra sin 9 sin fJ + 



2/2 2\ 

P (r -a ) 
rVA 



/A(, 



2r" 



(A4) 



(A5) 



7^ 



^ ^ (A(p^ -4r^) +p^(r^ - a^)) - 2ra sin sin /3 

- cos 9 cos asm 131 2a^A + (2(r^ + 

\ sin 9 ^ 

( ^ (r^ + a^) +a2s"'2« ^ 
sin a 2a V A + -^^ 



(A6) 



■ sin (3 



cos a sin (3 



r ( A + sin'' 



+ 2ra sin 9 sin /3 



(A7) 



These equations complete a first order set of 7 equations for the coordinates {t, r, 9, 0} and the parameters $, a and which 
is sufficient to determine photon's trajectory and is numerically very stable. A further advantage of this approach is that we 
have not yet made any use of the three conserved Killing quantities; the energy E, angular momentum J and Carter constant 
Q. In terms of spacetime position and our three photon momentum parameters, they are given by 



E = —(VA + a sin 9 sin 13) J = — (^asin^ 9 + sin 6* sin /3 (r^ + a^)) 
p p ^ \ / / 



<I>^ p^ ( 1 — cos^ a cos^ 



(A8) 



By evaluating these quantities explicitly at various points along a photon trajectory we can test that the integration is 
proceeding correctly. 

We apply the same basic idea to calculate the circular, timelike trajectories in the equatorial plane which model the flow 
of accreting matter, except now we can use a single parameter U , as follows: 



V — cosh Ut + sinh U (f> 

where the geodesic equations are satisfied by 
'r — a 



tanh U 



More commonly, these trajectories are written as 

1 



(A9) 



(AlO) 



(All) 
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and the constant of proportionality is obtained by imposing the condition v — 1. The two notations are related via 

n= g + tanh^^ (A12) 
r2 + + tanh U aVA 

In order to model the motion of a source orbiting with the disk matter, but at a general height hs — r sin above the equatorial 
plane, we can define a trajectory (no longer geodesic) as follows 

vs(xdt + nd^, n= ^. (ai3) 

a + (rsmOy''' 

effectively lifting a circular geodesic vertically out of the equatorial plane. The expression for U given by Equation lAlOl 
generalises to 

. , rr -a(r sin 6»)3/2 sinO 

acos^ + (r smO)-''^ -^A 

In order to measure the angle at which a photon leaves the source and hits the disk we must define a co-moving 3-frame 
attached to Vs (v being a special case where 6 — tv/2. Together with Equation ljA9^ . the vectors f, § and 0', where 

(j)' = sinh U(ji + cosh Ui, (A15) 

form a convenient orthonormal tetrad. This paper has been typeset from a T^jX/ I^T^riX file prepared by the author. 



